/*
 * FunctionListOld.cpp
 *
 *  Created on: Jun 27, 2011
 *      Author: fqhuy
 */
#include "FunctionsListOld.hpp"
#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <sstream>
//---------------------------------------------------------------------------
#define INS_FUNCS(T) \
		template void Deallocate3D(T ***arr3D, int l, int m); \
		template void Deallocate2D(T **arr2D, int rows); \
		template T** Allocate2D(int rows, int cols); \
		template T*** Allocate3D(int l, int m, int n); \
		template void solve_linear_eqn(T **A, T *B, int n, T *XX); \
		template std::string ToString3D(T ***array, const int& w, const int& h, const int& d,bool with_index); \
		template std::string ToString2D(T **array, const int& w, const int& h,bool with_index); \
		template int WriteArray3D(const char* file_name, T*** array, const int& w, const int& h, const int& d,bool with_index); \
		template int WriteArray2D(const char* file_name, T** array, const int& w, const int& h,bool with_index);

INS_FUNCS(float)
;
INS_FUNCS(double)
;
INS_FUNCS(int)
;
INS_FUNCS(long)
;
INS_FUNCS(char)
;

//allocate a 2D array
template<class T> T** Allocate2D(int rows, int cols) {
	T **arr2D;
	int i;

	arr2D = (T**) malloc(rows * sizeof(T*));
	for (i = 0; i < rows; i++) {
		arr2D[i] = (T*) malloc(cols * sizeof(T));
	}
	return arr2D;
}
template<class T> T*** Allocate3D(int l, int m, int n)
//void allocate3D(float ***arr3D,int l,int m,int n)
{
	T ***arr3D;
	int i, j, k;

	arr3D = (T***) malloc(l * sizeof(T **));

	for (i = 0; i < l; i++) {
		arr3D[i] = (T**) malloc(m * sizeof(T*));
		for (j = 0; j < m; j++) {
			arr3D[i][j] = (T*) malloc(n * sizeof(T));
		}
	}

	return arr3D;
}

template<class T>
void Deallocate3D(T ***arr3D, int l, int m) {
	int i, j;
	for (i = 0; i < l; i++) {
		for (j = 0; j < m; j++) {
			free(arr3D[i][j]);
		}
		free(arr3D[i]);
	}
	free(arr3D);
}

template<class T>
void Deallocate2D(T** arr2D, int rows) {
	int i;

	for (i = 0; i < rows; i++) {
		free(arr2D[i]);
	}
	free(arr2D);
}

template<class T>
void solve_linear_eqn(T **A, T *B, int n, T *XX) {
	T **AA, **L, **U;
	//AA = Fallocate2D(n, n + 1);
	//L = Fallocate2D(n, n);
	//U = Fallocate2D(n, n);

	AA = Allocate2D<T> (n, n + 1);
	L = Allocate2D<T> (n, n);
	U = Allocate2D<T> (n, n);

	// double XX[max_size]={0};

	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			AA[i][j] = A[i][j];

	for (int i = 0; i < n; i++)
		AA[i][n] = B[i];

	/*   for (int i = 0; i < n; i++)
	 for (int j = 0; j < n+1; j++)
	 ShowMessage(AA[i][j]);*/

	//----------------------  initialize_LU_matrices( )  --------------------//

	for (int count_1 = 0; count_1 < n; count_1++) {
		L[count_1][count_1] = 1;
		L[count_1][n] = AA[count_1][n];
	}

	for (int count_2 = 1; count_2 < n; count_2++)
		L[count_2][0] = (AA[count_2][0] / AA[0][0]);

	for (int count_3 = 0; count_3 < n; count_3++)
		U[0][count_3] = AA[0][count_3];

	//------------------  apply_LU_decomposition_method( )  -----------------//

	for (int count_1 = 0; count_1 < (n - 1); count_1++) {
		long double options[Max_Variables] = { 0 };

		int option_number = 0;

		for (int count_2 = (count_1 + 1); count_2 < n; count_2++) {
			for (int count_3 = 0; count_3 < n; count_3++)
				options[option_number] += (L[count_2][count_3]
						* U[count_3][(count_1 + 1)]);

			options[option_number] = (AA[count_2][(count_1 + 1)]
					- options[option_number]);
			option_number++;
		}

		int max_row_number = 0;

		for (int count_4 = 0; count_4 < option_number; count_4++) {
			if (fabs(options[max_row_number]) <= fabs(options[count_4])) {
				if (fabs(options[max_row_number]) == fabs(options[count_4])
						&& options[max_row_number] != options[count_4]) {
					if (options[max_row_number] < options[count_4])
						max_row_number = count_4;
				}

				else
					max_row_number = count_4;
			}
		}

		U[(count_1 + 1)][(count_1 + 1)] = options[max_row_number];

		long double temp = 0;

		for (int count_5 = (count_1 + 2); count_5 < n; count_5++) {
			for (int count_6 = 0; count_6 < n; count_6++)
				temp += (L[count_5][count_6] * U[count_6][(count_1 + 1)]);

			temp = AA[count_5][(count_1 + 1)] - temp;
			temp /= U[(count_1 + 1)][(count_1 + 1)];

			L[count_5][(count_1 + 1)] = temp;
			temp = 0;

			for (int count_7 = 0; count_7 < n; count_7++)
				temp += (L[(count_1 + 1)][count_7] * U[count_7][count_5]);

			temp = AA[(count_1 + 1)][count_5] - temp;
			U[(count_1 + 1)][count_5] = temp;
			temp = 0;
		}
	}

	//------------------------  generate_result( )  -------------------------//

	for (int count_1 = 0; count_1 < n; count_1++) {
		U[count_1][n] = L[count_1][n];

		for (int count_2 = 0; count_2 < count_1; count_2++)
			L[count_1][count_2] = (L[count_1][count_2] * U[count_2][n]);

		for (int count_3 = 0; count_3 < count_1; count_3++)
			U[count_1][n] -= L[count_1][count_3];

		U[count_1][n] /= L[count_1][count_1];
	}

	for (int count_4 = (n - 1); count_4 >= 0; count_4--) {
		XX[count_4] = U[count_4][n];

		for (int count_5 = (n - 1); count_5 > count_4; count_5--)
			U[count_4][count_5] = (U[count_4][count_5] * XX[count_5]);

		for (int count_6 = (n - 1); count_6 > count_4; count_6--)
			XX[count_4] -= U[count_4][count_6];

		XX[count_4] /= U[count_4][count_4];
	}
	Deallocate2D(AA, n);
	Deallocate2D(L, n);
	Deallocate2D(U, n);
}

template<class T> std::string ToString3D(T ***array, const int& w,
		const int& h, const int& d,bool with_index) {
	std::ostringstream oss(std::ostringstream::out);
	oss << std::fixed;

	for (int i = 0; i < w; i++) {
		for (int j = 0; j < h; j++) {
			for (int k = 0; k < d; k++) {
				if(with_index)
					oss<< i << " " << j << " " << k << " ";
				oss << array[i][j][k] << std::fixed << "\n";
			}
			//oss << std::endl;
		}
		//oss << std::endl;
	}

	return oss.str();
}

template<class T> std::string ToString2D(T **array, const int& w, const int& h,bool with_index) {
	std::ostringstream oss(std::ostringstream::out);
	oss << std::fixed;

	for (int j = 0; j < h; j++) {
		for (int k = 0; k < w; k++) {
			if(with_index)
				oss<< j << " " << k << " ";
			oss << array[j][k] << std::fixed << "\n";
		}
		//oss << std::endl;
	}
	return oss.str();
}

int WriteFile(const char* file_name, const char* str, const int& length) {
	FILE* fo;
	std::string fn = "output/";
	fn += file_name;

	fo = fopen(fn.c_str(), "wt");
	if (fo != NULL)
		for (int i = 0; i < length; i++)
			fputc(str[i], fo);
	else {
		printf("error occured while opening file %s\n", file_name);
		return -1;
	}
	fclose(fo);
	return 0;
}

template<class T> int WriteArray3D(const char* file_name, T*** array, const int& w, const int& h, const int& d,bool with_index){
	std::string data = ToString3D<T>(array,w,h,d,with_index);
	WriteFile(file_name,data.c_str(), data.length());
}

template<class T> int WriteArray2D(const char* file_name, T** array, const int& w, const int& h,bool with_index){
	std::string data = ToString2D<T>(array,w,h,with_index);
	WriteFile(file_name,data.c_str(), data.length());
}
